function [ stress, strain, rho ] = evaluateInteriorStress( F )
%EVALUATEINTERIORSTRESS Summary of this function goes here
%   Detailed explanation goes here

global poissonRatio shearModulus;

% F = evaluateInteriorDeformationGradient(x, y, domain2D, theta);

% F = -1 * F;

n = size(F, 3);

strain = zeros(2, 2, n);

stress = zeros(2, 2, n);

rho = zeros(n, 1);

II = [1, 0; 0, 1];

for i = 1 : n
    strain(:, :, i) = 0.5 * (F(:, :, i) + F(:, :, i)');
    stress(:, :, i) = shearModulus * (F(:, :, i) + F(:, :, i)') + 2 * shearModulus * poissonRatio / (1 - 2 * poissonRatio) * trace(F(:, :, i)) * II;
    rho(i) = tensorTensorContraction(strain(:,:,i), stress(:,:,i));
end

end

